home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / map_patch.pro < prev    next >
Text File  |  1997-07-08  |  9KB  |  224 lines
  1. ; $Id: map_patch.pro,v 1.12 1997/03/26 00:15:35 dave Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5.  
  6. function map_reduce_360, a    ;Reduce an angle to the range of +- 180.
  7. b = a mod 360.
  8. neg = where(b lt 0, count)
  9. if count gt 0 then b[neg] = b[neg]+360.
  10. return, b
  11. end
  12.  
  13.  
  14. FUNCTION  map_patch, Image_Orig, Lons, Lats, $
  15.         XSTART = xstart, YSTART = ystart, $
  16.         XSIZE = xsize, YSIZE = ysize, $
  17.         LON0 = lon0, LON1 = lon1, $
  18.         LAT0 = lat0, LAT1 = lat1, $
  19.         MISSING = missing, MAX_VALUE = max_value, $
  20.         TRIANGULATE=triangulate, DEBUG = debug
  21. ;+
  22. ;NAME:
  23. ;     map_patch
  24. ;PURPOSE:
  25. ;    This function linearly interpolates a data sampled in
  26. ;    latitude/longitude into device space.  Data values may be
  27. ;    either rectangularly or irregularly gridded over the globe.
  28. ;Category:
  29. ;        Mapping
  30. ;Calling Sequence:
  31. ;        result = map_patch(Image_Orig [, Lons] [, Lats])
  32. ;INPUT:
  33. ;      Image_Orig- A array containing the data to be overlayed on map.
  34. ;        It may be either 1D or 2D.  If 2D, it has Nx columns,
  35. ;        and Ny rows.  The cell connectivity must be
  36. ;        rectangular, unless the TRIANGULATE keyword is specified.
  37. ;        If Image_Orig is 1D, then Lats and Lons must contain
  38. ;        the same number of points.
  39. ;    Lons-    A vector or 2D array containing the longitude of each
  40. ;        data point or column.  If Lons is 1D, lon(image_orig(i,j)) =
  41. ;        Lons(i); if Lons is 2D, lon(image_orig(i,j)) = Lons(i,j).
  42. ;        This optional parameter may be omitted if the
  43. ;        longitudes are equally spaced and are
  44. ;        specified with the LON0 and LON1 keywords.
  45. ;    Lats-    A vector or 2D array containing the latitude of each
  46. ;        data point or row  If Lats is 1D, lat(image_orig(i,j))
  47. ;        = Lats(j); if Lats is 2D, lat(image_orig(i,j)) =
  48. ;        Lat(i,j). This optional parameter may be omitted
  49. ;        if the latitudes are equally spaced and are specified
  50. ;        with the LAT0 and LAT1 keywords.
  51. ;KEYWORDS:
  52. ;    LAT0-    The latitude of the first column of data.  Default=-90.
  53. ;    LAT1-    The latitude of the last column of data.  Default=+90.
  54. ;    LON0-    The longitude of the first row of data.  Default=-180.
  55. ;    LON1-    The longitude of the last row of data.  Default=180-360/Ny.
  56. ;    MISSING = value to set areas outside the valid map coordinates.
  57. ;        If omitted, areas outside the map are set to 255 (white) if
  58. ;        the current graphics device is PostScript, or 0 otherwise.
  59. ;    MAX_VALUE = values in Image_Orig greater than MAX_VALUE
  60. ;        are considered missing.  Pixels in the output image
  61. ;        that depend upon missing pixels will be set to MISSING.
  62. ;    TRIANGULATE = if set, the points are Delauny triangulated on
  63. ;        the sphere using the TRIANGULATE procedure to determine 
  64. ;         connectivity.    Otherwise, rectangular connectivity is assumed.
  65. ; Optional Output Keywords:
  66. ;    xstart --- the  x coordinate where the left edge of the image
  67. ;        should be placed on the screen.
  68. ;    ystart --- the y coordinate where th bottom edge of the image
  69. ;        should be placed on the screen.
  70. ;    xsize ---  returns the width of the resulting image expressed in
  71. ;        graphic coordinate units.  If current graphics device has
  72. ;        scalable pixels,  the value of XSIZE and YSIZE should
  73. ;        be passed to the TV procedure.
  74. ;    ysize ---  returns the pixel height of the resulting image, if the
  75. ;        current graphics device has scalable pixels. 
  76. ;
  77. ; Restrictions:
  78. ;    This could be quicker.
  79. ; Output:
  80. ;      The interpolated function/warped image is returned.
  81. ;
  82. ; Procedure:  
  83. ;    An object space algorithm is used, so the time required 
  84. ;    is roughly proportional to the size, in elements, of Image_Orig.
  85. ;    Computations are performed in floating point.
  86. ;    For rectangular grids, each rectangular cell of the original
  87. ;    image is divided by a diagonal, into two triangles.  If
  88. ;    TRIANGULATE is set, indicating irregular gridding, the cells are
  89. ;    triangulated.  The trianges are then converted from lat/lon to  
  90. ;    image coordinates and then interpolated into
  91. ;    the image array using TRIGRID.
  92. ;
  93. ;MODIFICATION HISTORY:
  94. ;    DMS of RSI, July, 1994.        Written.
  95. ;    DMS, Nov, 1996.    Rewritten and adapted for new maps, rev 2.
  96. ;-
  97.  
  98. ON_ERROR,2
  99.  
  100. if (!x.type NE 2) and (!x.type ne 3) THEN  $        ;Need Mapping Coordinates
  101.   message, "Current window must have map coordinates"
  102.  
  103. s = size(Image_Orig)
  104. Nx = s[1]                       ; # of columns
  105. if s[0] eq 2 then Ny = s[2] else Ny = 1 ; # of rows
  106. n = N_elements(Image_orig)
  107.  
  108. if n_elements(lons) eq 0 then begin ;Make longitudes?
  109.     if n_elements(lon0) le 0 then lon0 = -180.
  110.     if n_elements(lon1) le 0 then lon1 = 180. ;As documented
  111. ;    if n_elements(lon1) le 0 then lon1 = lon0 - 360./nx + 360. ;as it was
  112.     dx = lon1-lon0
  113.     if dx le 0 then dx = dx + 360.
  114.     lons = findgen(nx) * (dx/(nx-1.)) + lon0
  115. endif
  116.  
  117. if n_elements(lats) eq 0 then begin ;Make lats?
  118.     if n_elements(lat0) le 0 then lat0 = -90.
  119.     if n_elements(lat1) le 0 then lat1 = 90.
  120.     lats = findgen(ny) * ((lat1-lat0)/(ny-1.)) + lat0
  121. endif
  122.  
  123. if n_elements(lats) ne ny and n_elements(lats) ne n then $
  124.   message, "Lats has incorrect size"
  125. if n_elements(lons) ne nx and n_elements(lons) ne n then $
  126.   message, "Lons has incorrect size"
  127.  
  128. lonmin = min(lons, MAX=lonmax)
  129. latmin = min(lats, MAX=latmax)
  130. wrap = abs(map_reduce_360(lonmax + 360./nx-lonmin)) lt 1e-4
  131. if wrap eq 0 then lnlim = [lonmin, lonmax] else lnlim = [-180., 180.]
  132. ltlim = [latmin, latmax]
  133. scale = !d.flags and 1        ;TRUE if device has scalable pixels
  134.  
  135. IF scale THEN BEGIN        ;Fudge for postscript?
  136.     scalef = 0.02        ;PostScript scale factor
  137.     scale_orig = [!x.s, !y.s]
  138.     !x.s = !x.s * scalef
  139.     !y.s = !y.s * scalef
  140. ENDIF ELSE scalef = 1
  141.  
  142. xw = scalef * !x.window         ;Map extent in normalized coords
  143. yw = scalef * !y.window
  144.  
  145. xw1 = !map.uv_box[[0,2]] * !x.s[1] + !x.s[0]
  146. yw1 = !map.uv_box[[1,3]] * !y.s[1] + !y.s[0]
  147.  
  148. ;Screen extent of our image in device coords
  149. screen_x = long([ xw[0] > xw1[0], xw[1] < xw1[1] ] * !d.x_size)
  150. screen_y = long([ yw[0] > yw1[0], yw[1] < yw1[1] ] * !d.y_size)
  151.  
  152. if n_elements(max_value) eq 0 then max_value = max(Image_orig)
  153. if n_elements(missing) le 0 then missing = 0
  154.  
  155. rect = [screen_x[0], screen_y[0], screen_x[1], screen_y[1]]
  156.  
  157. if keyword_set(triangulate) then begin ;Make our own triangulation?
  158.     if n_elements(lons) eq n_elements(image_orig) then x = lons $
  159.     else if n_elements(lons) eq nx then x = lons # replicate(1., ny) $
  160.     else message, 'Dimension of LONS incompatible'
  161.     
  162.     if n_elements(lats) eq n_elements(image_orig) then y = lats $
  163.     else if n_elements(lats) eq ny then y =  replicate(1.,nx) # lats $
  164.     else message, 'Dimension of LATS incompatible'
  165.     
  166.     temp = image_orig           ;Copy cause TRIANGULATE, /SPHERE reorders data
  167.     triangulate, x, y, triangles, Sphere=s, /DEGREES, FVALUE=temp
  168.     x = trigrid(x, y, temporary(temp), triangles, [1., 1.], rect, $
  169.                 MAX_VALUE=max_value, MISSING = missing, $
  170.                 MAP=[!x.s * !d.x_size, !y.s * !d.y_size])
  171.     
  172. endif else begin        ;Make a regular triangular grid
  173.     
  174.                                 ; Make 2D lat/lon arrays
  175.     if n_elements(lons) eq nx then x = lons # replicate(1., ny)
  176.     if n_elements(lats) eq ny then y = replicate(1.,nx) # lats
  177.     nr = nx - 1 + wrap        ;If wrap, connect last row to first
  178.     do_north = abs(latmax-90.) lt 1.0e-3
  179.     do_south = abs(latmin+90.) lt 1.0e-3
  180. ;        # of rows of triangles, 2 / row, except at poles.
  181.     nr1 = 2 * (ny-1) - do_north - do_south
  182.     i = (nx-1) * 6              ;Vertices / row
  183.     t = lonarr(nr * 6, /NOZERO) ;Typical row of rects
  184. ; Make our triangles in COUNTERCLOCKWISE order so they'll be properly clipped
  185. ; and split by the map routines.
  186.     for i=0, nr-1 do t[i*6] = ([i,i+1,  i, i,i+1,i+1] mod nx) + $
  187.       [0, nx, nx, 0,  0, nx]
  188.     triangles = lonarr(nr1 * nr * 3, /NOZERO)
  189.     j = 0L
  190.     y0 = 0
  191.     if do_south then begin
  192.         for i=0, nr-1 do triangles[i*3] = ([i,i+1,i] mod nx) + [0, nx, nx]
  193.         j = 3*nr
  194.         y0 = 1
  195.     endif
  196.     if do_north then y1 = ny-3 else y1 = ny-2 ;Last row of rects
  197.     for iy= y0, y1 do begin
  198.         triangles[j] = t + iy * nx
  199.         j = j + 6*nr
  200.     endfor
  201.     if do_north then $          ;Top row
  202.       for i=0, nr-1 do triangles[j+i*3] = ([i, i+1, i] mod nx) + $
  203.       [0,0, nx] + nx*(ny-2)
  204.     
  205. ; Interpolate from the triangular grid, in lat/lon onto an image bitmap.
  206.     x = trigrid(x, y, image_orig, triangles, [1., 1.], rect, $
  207.                 MAX_VALUE=max_value, MISSING = missing, $
  208.                 MAP=[!x.s * !d.x_size, !y.s * !d.y_size])
  209. endelse                         ;triangular grid
  210.  
  211.  
  212. xsize = ceil((rect[2] - rect[0])/scalef)
  213. ysize = ceil((rect[3] - rect[1])/scalef)
  214.  
  215. xstart = long(rect[0] / scalef)
  216. ystart = long(rect[1] / scalef)
  217.  
  218. if scale then begin        ;Restore scale
  219.     !x.s = scale_orig[0:1] & !y.s = scale_orig[2:3]
  220.     endif
  221.  
  222. return, x
  223. end
  224.